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We present two new adaptive quadrature routines. Both routines differ from previously published 
algorithms in many aspects, most significantly in how they represent the integrand, how they 
treat non-numerical values of the integrand, how they deal with improper divergent integrals and 
how they estimate the integration error. The main focus of these improvements is to increase the 
reliability of the algorithms without significantly impacting their efficiency. Both algorithms are 
implemented in Matlab and tested using both the "families" suggested by Lyness and Kaganove 
and the battery test used by Gander and Gautschi and Kahaner. They are shown to be more 
reliable, albeit in some cases less efficient, than other commonly-used adaptive integrators. 

Categories and Subject Descriptors; F.2.1 [Numerical Analysis]: Numerical Algorithms and 
Problems — Computations on polynomials; G.1.4 [Numerical Analysis]: Quadrature and Nu- 
merical Differentiation — Adaptive and iterative quadrature 

General Terms: Algorithms, Performance, Reliability 

Additional Key Words and Phrases: Adaptive quadrature, interpolation, orthogonal polynomials, 
error estimation 



1. INTRODUCTION 

Since the publication of the first adaptive quadrature algorithms almost 50 years 
ago, much has been done and even more has been written on the subject^: In 
1962 Kuncir [Kuncir 1962] kickcd-off the ficld^ with his adaptive Simpson's rule 
integrator, which uses - as the name suggests - Simpson's rule to approximate 
the integral, bisecting recursively until the difference between the approximation in 
one interval and that of its two sub-intervals is below the required tolerance. This 



^A recent review by the author [Gonnet 2009a] limited to error estimation listed 21 distinct 
published algorithms and references to more than 50 publications directly related to adaptive 
quadrature. 

^Davis and Rabinowitz [Davis and Rabinowitz 1984] reference, as the first adaptive quadrature 
routines, the works of Villars [Villars 1956], Henriksson [Henriksson 1961] and Kuncir [Kun- 
cir 1962]. Henriks.son'H algorithm, which appeared in the first issue of BIT, is an ALGOL- 
implementation of the algorithm described by Villars, which is itself an extension of an algorithm 
by Morrin [Morrin 1955]. These algorithms, however, are more reminiscent of ODE integrators, 
which is why we will not consider them to be "genuine" adaptive quadrature routines. 
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approach, as simple as it may seem, still lives on, with some minor modifications, as 
the default integrator quad in MATLAB (added in [The Mathworks 2003]). quad is 
itself a modification of Gander and Gautschi's 2001 adaptsim [Gander and Gautschi 
2001] which is itself a modification of Lyness' 1970 SQUANK [Lyness 1970], which is 
itself almost identical to Kuncir's original algorithm^. 

At the time of Kuncir's publication, McKeeman [McKeeman 1962; McKeeman 
and Tesler 1963; McKeeman 1963] and later Forsythe et al. [Forsythe et al. 1977] ex- 
tended this approach to use higher-degree Newton-Cotes rules and/or sub-division 
into more than two sub-intervals. 

In 1971, de Boor [de Boor 1971] introduced the concept of "double adaptivity", 
constructing a Romberg T-table [Bauer et al. 1963] within each sub-interval to 
approximate the integrand and deciding, at every step, whether to extend the T- 
table by another row {i.e. increase the order of the quadrature) or to subdivide the 
interval. This decision was made by using the convergence rates of the columns 
of the T-table to guess the integrand's behavior, i.e. "well behaved", singular, 
discontinuous or noisy, and apply an adequate strategy for that behavior. 

The emergence of more powerful computers and better algorithms {e.g. Golub 
and Welsch [1969] and Gentleman [1972]) for the construction of more complex 
quadrature rules quickly led to the wider use of Clenshaw-Curtis quadrature rules 
[Clenshaw and Curtis 1960], as used by O'Hara and Smith [1969] and Oliver [1972], 
and later the use of Gauss quadrature rules and their Kronrod extensions [Kronrod 
1965], first used by Piessens [1973] and Patterson [1973] independently. 

Other interesting and/or noteworthy advances in the field are: 

— The introduction of stratified or recursively monotone stable (RMS) quadrature 
rules [Sugiura and Sakurai 1989; Favati et al. 1991; Laurie 1992], filling the gap 
between low-degree (due to low numerical stability at higher degrees) Newton- 
Cotes rules, the nodes of which nodes are re-usable over several recursion levels, 
and high-degree (due to better numerical stability) yet non-reusable Clenshaw- 
Curtis or Gauss rules, thus providing some extra efficiency, 

— The use of non-linear extrapolation when computing the integral or the error esti- 
mate [Rowland and Varol 1972; Venter and Laurie 2002; Laurie 1983; de Doncker 
1978], as is done in the highly successful QAGS subroutine in the QUADPACK in- 
tegration library, 

— The use of higher-order coefficients relative to some base to compute the error 
estimate [O'Hara and Smith 1968; Oliver 1972; Berntsen and Espclid 1991]. 

A number of authors have published comparisons of these and many other adap- 
tive quadrature routines [Casaletto et al. 1969; Hillstrom 1970; Kahaner 1971; Mal- 
colm and Simpson 1975; Robinson 1979; Krommcr and Ubcrhubcr 1998; Gonnet 
2009a] as well as methodologies to compare different routines [Lyness and Kaganove 
1977]. 

As already noted in Rice [1975], despite all their differences, most adaptive 
quadrature algorithms follow the general scheme, as in Algorithm 1. First, an 



^Lyness himself formulates his algorithm as an extension to McKeeman's Adaptive Integrator [Mc- 
Keeman 1962], yet the resulting algorithm is much more similar to Kuncir's, which was published 
merely a few months before McKeeman's. 
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estimate of the integral in the interval [a,b] is computed (Line 1). An error esti- 
mate of the integral is then computed (in this example, an absolute error estimate 
is approximated, Line 2). If this estimate is smaller than the required tolerance 
(Line 3), then the estimate is returned (Line 4). Otherwise, the interval is bisected 
and the algorithm is called on both halves (Line 6) using a modified local tolerance 
t'. 



Algorithm 1 int (/, a, b, r) 



^~ - /a /(a;) da; 
if £ < r then 
return Q 
else 

return int(/, a, (a + b) /2, r') + int(/, (a + 6)/2, 6, t') 

{call the integrator recursively on both sub-intervals) 

end if 



{approximate the integral in [a,b]) 
{approximate the integration error) 

{return the current estimate) 



Algorithm 2 int (/, a, b, r) 



1 

2 

3 
4 
5 
6 
7; 
8 

9 
10 
11 
12 



Qo « /„ f{x) dx 

Q-j;;f{x)dx 

{[a,b, Qo,so]} 



So 
H 

while ^ 



e, > r do 



k -ir- arg maxfe Sk 

H H\{[ak,bk,Qk-,£k\} 
mJ^{ak + bk)/2 

Qieft « /™ fix) 

Qieft - Q f{x) dx 



eieft 



•bright « !^ f{x) dx 



^right 



{approximate the integral in [a,b]) 
{approximate the integration error) 
{initialize the heap with the first interval) 

{pop the interval with the largest error) 



{compute the integral on the left) 
{compute the error on the left) 
{compute the integral on the right) 
Qright — J^ji f{x) dx {compute the error on the right) 

HU {[afc,m, (3left,£leit], ["l, ^'fe, <5 right, £right]} 

{push the new intervals back on the heap) 



13; end while 
14: return J2 



Qi 



{return the sum of the integrals in the intervals) 



Not all adaptive quadrature algorithms are recursive (locally adaptive): many 
algorithms, such as those in QUADPACK, maintain a heap of intervals and bisect 
the interval with the largest local error estimate and return the new sub-intervals 
to the heap until the sum of the local errors is below the required tolerance (Algo- 
rithm 2, globally adaptive). This approach, although more memory-intensive, has 
several advantages over the recursive approach, such as better control over the error 
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estimate and the ability to restart or refine an initial approximation [Malcolm and 
Simpson 1975; Rice 1975]. 

However, despite all these advances in numerical quadrature in general and adap- 
tive quadrature specifically, the results of these methods must often be treated with 
caution, as failures are common even for relatively simple integrands. In this paper 
we will present two new adaptive quadrature algorithms which attempt to address 
this lack of reliability. The algorithms follow the general scheme in Algorithm 2, 
yet with significant differences to previous methods regarding how the integrand 
is represented (Section 2), how the integration error is estimated (Section 3) and 
how singularities (Section 4) and divergent integrals (Section 5) are treated. The 
algorithm itself is presented in Section 6 and in Section 7 it is validated against 
other popular algorithms. These results are then discussed in Section 8. 

2. FUNCTION REPRESENTATION 

In most quadrature algorithms, the integrand is not represented internally except 
through different approximations of its integral. We denote such approximations 

as 



where n is the degree^ of the quadrature rule and m its multiplicity. 

The quadrature rule Q„ itself is computed as the weighted sum of the integrand 
evaluated at a pre-determined set of nodes^ Xj € [— 1, 1], i = . . . n: 



The evaluation of one or more such quadrature rules is usually the only informa- 
tion considered regarding the integrand. 

Some authors [Gallaher 1967; Ninomiya 1980] use additional nodes to numerically 
approximate the higher derivative directly using divided differences, thus supplying 
additional information on the integrand f{x). In a similar vein, O'Hara and Smith 
[1969], Oliver [1972] and Berntsen and Espelid [1991] compute some of the higher- 
order coefficients of the function relative to some orthogonal base, thus further 
characterizing the integrand. 

In all of these cases, however, the characterization of the integrand is not com- 
plete and in most cases only implicit. In the following, we will attempt to better 
characterize the integrand. 

Before doing so, we note that for every interpolatory quadrature rule, we are in 



*In the following, we will use the term "degree" to specify the algebraic degree of precision of a 
quadrature rule, which is the highest degree for which all polynomials of that degree will always 
be integrated exa<;tly by the rule. 

^In the following we assume, for notational simplicity, that the number of nodes is the degree of 
the rule plus one. Although most quadrature rules, e.g. interpolatory quadrature rules with an 
odd number of symmetric nodes or Gauss quadratures and their Kronrod extensions, need less 
than n + 1 nodes for degree n, this is a general upper bound for interpolatory quadrature rules. 
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fact computing a interpolating polynomial gn{x) of degree n such that 

gnixi) = f{xi), i = 0...n 
and evaluating the integral thereof 



Q„[a,b] = / gn{x)dx. 

J a 

This ciquivalence is easily demonstrated, as is done in many textbooks in numeri- 
cal analysis ([Stiefel 1961; Rutishauser 1976; Gautschi 1997; Schwarz 1997; Ralston 
and Rabinowitz 1978] to name a few)^. 

Since any polynomial interpolation of dcigree n over n + 1 distinct points is 
uniquely determined, it doesn't matter how we choose to represent gn{x) - its 
integral will always be identical to the result of the interpolatory quadrature rule 
Qn[a, b] over the same nodes. 

In the following, we will represent as a linear combination of orthogonal 

polynomials: 

n 

9nix) =^CiPi{x) (2) 

where the pi{x), i = . . .n are polynomials of degree i which are orthonormal with 
respect to some inner product 

j^k, 
j = k. 

We will use the coefficients c = (cq, ci, . . . , c„)^ from Equation (2) as our represen- 
tation of 5,1 (x)- 

For notational simplicity, we will assume that the integrand has been transformed 
from the interval [a,b] to the interval [—1,1]- The polynomial gn{x) interpolates 
the integrand f{x) at the nodes Xi G [—1, 1]: 

gn{xi) = f{xi), i = 0...n. 

Given the function values f = (/(.xo). /(xi). . . . , /(.Tn))''' at the nodes Xi, i = 
. . . n, we can compute the coefficients by solving the linear system of equations 

Pc = f (3) 

where the matrix P with = Pj{xi) on the left-hand side is a Vandermonde- 
like matrix. The naive solution using Gaussian elimination is somewhat costly and 
may be unstable [Gautschi 1975]. However, several algorithms exist to solve this 
problem stably in 0{n^) operations for orthogonal polynomials satisfying a three- 
term recurrence relation [Bjorck and Pereyra 1970; Higham 1988; 1990; Gonnet 
2009b]. 



{pj,Pk) = 1 5 



®If we consider the Lagrange interpolation g„{x) of the integrand and integrate it, we obtain 

Zb rb " ^ rb " 

g„{x)dx= / '^ei{x)f{xi)dx = '^f{xi) li{x) dx = y^J(xi)wi 
i=0 i=0 1=0 

where the li {x) are the Lagrange polynomials and the Wi are the weights of the resulting quadra- 
ture rule. 
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In the following, we will use the orthonormal Legendre polynomials, which are 
orthogonal with respect to the inner product 



,Pk) = j 



Pj{x)pk{x) dx. 



(4) 



We will evaluate and interpolate the integrand at the Chebyshev nodes 



COS 



n 



0...n. 



These nodes are chosen over the Gauss quadrature nodes or equidistant nodes due 
to their stability [Trcfcthcn 2008] , because the nodes can be rc-uscd whcin increasing 
the degree of the rule [Oliver 1972] and because they include the interval boundaries. 

The resulting Vandermonde-like matrix has a condition number (P) e 0(n3/2) 
which remains < 1000 for n < 100 and is thus tractable even for moderate n 
[Gonnet 2009a]. 

The resulting representation of gn{x) (Equation (2)) has some interesting prop- 
erties. First of all, it is simple to evaluate the integral of gn{x) using 



/I ^1 n ri .1 

gn{x)dx= I V'ciPi(a;)da; = V'cj / Pi{x)dx = u)^ 



(5) 



where the weights can be pre-computed and applied much in the same way as the 
weights of a quadrature rule. Note that for the normalized Legendre polynomials 
used herein, = (l/-\/2, 0, . . . , 0). 

We can also evaluate the L2-norm of gn{x) quite efHciently using Parseval's the- 
orem 



gl{x) dx 



1/2 



.1=0 



1/2 



In the following, we will use || • || to denote the 2-norm for vectors. 

A final useful feature is that, given the coefficients of gn{x) on [—1,1], we can 
construct upper-triangular matrices 

= J_^Pi{x)Pj (^) da;, tI;} = J^^Pi{x)pj (^) dx, i = 0...n,j>i 



such that 



cW=tWc and cW=tWc 



(6) 



are the coefficients of gn{x) on the left and right sub-intervals [—1,0] and [0, 1] re- 
spectively. These matrices depend only on the polynomials Pi{x) and can therefore 
be pre-computed for any set of nodes such that'^ 



1 



xe [-1,1]. 



(7) 



''This can be shown by representing the polynomials Pi({x — l)/2) of degree i as a linear combi- 
nation of the polynomials Pj{x), j = . . .i where the coefficients are computed using the inner 
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The coefficients c'^^ and c'^''^ can be useful if, after bisecting an interval, we want to 
re-use, inside one of the sub-intervals, the interpolation computed over the entire 
original interval. 

3. ERROR ESTIMATION 

This section contains a summary of the more important results of [Gonnct 2009a]. 
For a more complete discussion and testing of the error estimators presented herein, 
we refer to that publication. 

Although they differ in their specific implementations, what all these error esti- 
mates have in common is that they try to approximate the quantity 



Q<rHa,b]- f f{x)dx 

J a 



(8) 



using only two or more approximations of the integral or of its coefficients relative 
to some base. In these estimates, problems may occur when the difference between 



two estimates (5„ [a, 6] and (5„ ^ [a, 6] or Q„j [a, 6] and Qnal^i^li or the magni- 
tude of the computed coefficients is accidentally small^, i.e. the approximations 
used to compute the error estimate are too imprecise, resulting in a false small er- 
ror estimate. This is often the case near singularities and discontinuities where the 
assumptions on which the error estimate is based, e.g. continuity and/or smooth- 
ness, do not hold. If we re-construct the underlying interpolatory polynomials for 
the pair of quadrature rules used, we see that the interpolations differ significantly 
{e.g. see Figure 1). This difference is a good indicator for whether the integrand is 
correctly represented or not. 

It is for this reason that instead of trying to compute the error as in Equation (8), 
we will use the L2 norm^ of the difference between the integrand f{x) and the 



Pj{x)Pi 



product in Equation (4): 

We can then re-insert this expression into Equation (7) 

'x-1 



i r /■! 



9n^ (x) = gn 



i=0 



dx 



x~l 



j=0 



3>» ■ 



i=0 j=0 



which, swapping the indices i and j and re-arranging the sums on the right-hand side can be 
re-written as 



j2cfp,ix) = j:p,ix)j2cK 



using which the vector of coefficients c(^) can be computed as in Equation (6). 

®This term was first used by O'Hara and Smith [1968] to describe this problem. 

®Note that if instead of using the orthonormal Legendre polynomials we were to use polynomials 

orthogonal with respect to any specific measure w{x), we would compute the L2-norm with respect 

to that measure: 

- rl -1 1/2 

£= J w{x){gn{x) ~ fix))"^ dx 

Fortunately enough, for any measure 'w{x), the following derivations apply without modification. 
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Fig. 1. A discontinuous function (solid) and the integrand interpolation using the Gauss-Lobatto 
rule Q5[— 1,1] (dashed) and the Gauss-Kronrod rule (59[— 1, 1] (dotted) from Matlab's quadl in- 
tegrator. Note that although both quadratures return the same result, the interpolated functions 
differ significantly. 



interpolant gn{x) 



{9n{x)-f{x)) dx 



1/2 



(9) 



This is an approximation of the integration error Equation (8). The error estimate 
will only be zero if the interpolated integrand matches the integrand on the entire 
interval 

9n{x) = /(x), X e [-1, 1]. 

In such a case, the integral will also be computed exactly. The error Equation (9) 
is therefore, assuming we can evaluate it reliably, not susceptible to "accidentally 
small" values. 

Since we do not know f{x) explicitly, i.e. we can only sample f{x) in a point- wise 
fashion, we cannot evaluate the right-hand side of Equation (9) exactly. In a first, 
naive approach, we could compute two interpolations gn^ (x) and gn^ (x) of different 
degree where ni < 712. If we assume, as is done for error estimators using quadrature 
rules of differing degree, that gn2 (x) is a sufRciently precise approximation of the 
integrand 

gl^Jix)- fix), xe[-i,l] 

then we can approximate the error of the interpolation gl^^ (x) as 



fix)-g^^Hx) 



-.1/2 



da; 



{x)-gi!^Hx)y 



dx 



1/2 



= ||c(2)-c(i)| 
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where c*^^^ and c*^^^ are the vectors of the coefficients of gni{x) and Qn-^ix) respec- 
tively and = where i> ni. Our first, naive error estimate is hence 



^na 



,(1) 



,(2) 



(10) 



This error estimate, however, only applies to the lowcr-dcgrcc estimate g^}{x). 
Yet if we are going to compute the higher-degree estimate gn^ix), it would be 
preferable to have an error estimate for that approximation. 

We can, taking a different approach, use the interpolation error 



\9n{x) - f{x)\ = 



(n + l)! ""(^^ 



, xe[-i,i] 



(11) 



where G [—1,1] depends on the value of x and where 7r„(a;) = n"=o(*^ — Xi) is 
the Newton basis polynomial over the nodes of the interpolation gn{x). Taking the 
i2-iiorm on both sides of Equation (11) we obtain 



£ = 



f{x)f dx 



1/2 



(n + l)! 



-,1/2 



nl{x) dx 



Since Tr^{x) is, by definition, positive for any x, we can apply the mean value 
theorem of integration and extract the derivative resulting in 



j ^{9n{x) - f{x)f dx 



1/2 



(n+l)! 



j ^-^lix) 



dx 



1/2 



, ^e[-i,i]. 

(12) 

Given two interpolations gn\x) and gn\x) over a non-identical set of nodes, we 
can compute the interpolation errors 



ii:\x)-f{x) 



(n+l)! 



<*'{x) 



e[-i,i], *G{1,2} (13) 



where iTn^ {x) and TVn^ (x) are the Newton basis polynomials over the nodes of gn^ {x) 
and gn\x) respectively. If wc assume that /^"^^•'(.x) is constant for x G [—1, 1] (a 
stricter version of the "sufficiently smooth" assumption for the purpose of deriving 
this error estimate) and take the i2-iiorm of the difference between both errors, we 
obtain 



'^(g^:\x)-g(:\x))' dx 



,1/2 



J(n+1)(^) 



(n+l)! 



'^(7rW(x)-7r(^)(x))^ 



da; 



1/2 



(14) 



If we represent the interpolations gn^ [x) and gn^ {x) by their coefficients c^^^ and 
c^^) respectively, then we can write the left-hand side of Equation (14) as 



,(1) 



,(2) 



J(n+1)(^) 



(n + l)! 



y'^(7rW(x)-42)(x))'dx 



Similarly, if we represent the Newton basis polynomials 'Kn\x) and iTn^x) by their 
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coefficients b*^^) and b'^' respectively 



r(l) 



(x) 



n+1 

i=0 



n+1 



we can isolate the fraction on the right hand side of Equation (15) 



(n + 1)! 



.(1) 



b(l) -b(2) 



(16) 



(17) 



Inserting this expression into the original error estimate (Equation (12)) for the 
interpolation Qn^ (x) we then obtain 



dx 



= llbWl 



b(i) -b(2) 



=: £(1) 



Hence, using two interpolations of the same degree, we obtain the more refined 
error estimate 



Sref : = 



llcW -c(2)| 



llbWl 



(18) 



||b(l) -b(2)| 

for the interpolation gn\x). 

Note that if the nodes of the interpolations gn^ (x) and gli^ (x) are fixed, we can 
pre-compute the scaling Hb^^^ ||/||b'-^' — b^^^H- 

Instead of explicitly computing two different interpolations over an interval [a, b] 
to construct the error estimate, we can rc-\isc the interpolation from the previoiis 
level of recursion after bisection, the coefficients c°''' of which can be computed 
using Equation (6). Likewise, we can compute the coefficients b"''' of the Newton 
basis polynomial over the nodes of the previous level in the same way. However, 
since b°''' and b are not in the same interval, we have to scale the coefficients of b"'*^ 
by 2"+^ such that Equation (13) holds. In this way, we only need to compute and 
store a single matrix and vector b for a single stencil of interpolation nodes. 

As with the previous error estimators, we have also made an assumption of 
smoothness regarding the integrand by assuming that /("+^^(a;) is constant for 
x e [—1,1] to construct Equation (18). We can't verify this directly, but we can 



verify if our computed | - 



(«+i) 



(ill 



(Equation (17)) actually satisfies Equation (13) 



for the nodes of the first interpolation by testing if 



fiXi) 



(n + 1)! 



r(2) 



(Xi) 



(19) 



is satisfied for alH = . . . n, where the Xj are the nodes of the interpolation (x). 
The value i^i > 1 is an arbitrary relaxation parameter (for the tests in Section 7 we 
use i^i = 1.1). If this condition is violated for any of the Xi, then we use the naive 
error estimate in Equation (10). 

4. SINGULARITIES AND UNDEFINED VALUES 

Since most adaptive quadrature algorithms are designed for general-purpose use, 
they will often be confronted with integrands containing singularities or undefined 
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function values. These can cause problems on two levels: 

— The quadrature rule has to be evaluated with a non-numerical value such as a 
NaN or ±lnf, 

— The integrand may not be as smooth and continuous as the algorithm might 

assume. 

Such problems arise when integrating functions such as 

rh 

dx, a <0 



J 

Jo 



which have a singularity at a; = 0, or when computing seemingly innocuous integrals 
such as 

''^ sinx , 
da; 

X 

for which the integrand is undefined at x = 0, yet has a well-defined limit 

sin X 
lim = 1. 

x^O X 

In both cases problems could be avoided by either shifting the integration domain 
slightly or by modifying the integrand such as to catch the undefined cases and 
return a correct numerical result. This would, however, require some prior reflection 
and intervention by the user, which would defeat the purpose of a general-purpose 
quadrature algorithm. 

Most algorithms deal with singularities by ignoring them, setting the offending 
value of the integrand to [Davis and Rabinowitz 1984, Section 2.12.7]. Another 
approach, taken by quad and quadl in Matlab, is to shift the edges of the domain 
by Emach if non-mimcrical value is encoimtcred there and to abort with a warning 
if a non-numerical value is encountered elsewhere in the interval. Since singularities 
may exist explicitly at the boundaries {e.g. integration of x", a < in the range 
[0, /i]), the explicit treatment of the boimdaries is needed, whereas for arbitrary 
singularities within the interval, the probability of hitting them exactly is somewhat 
small. 

QUADPACK's QAG and QAGS algorithms take a similar approach: since the nodes 
of the Gauss and Gauss-Lobatto quadrature rules used therein do not include the 
interval boundaries, non-numerical values at the interval boundaries will be implic- 
itly avoided. If the algorithm has the misfortune of encountering such a value inside 
the interval, it aborts. 

Our approach to treating singularities will be somewhat different: instead of 
setting non-numerical values to 0, we will simply remove that node from our in- 
terpolation of the integrand. This can be done rather efficiently by computing the 
interpolation as shown before using a function value of f{xj) = for the offending 
jth node and then down-dating (as opposed to up-dating) the interpolation, i.e. re- 
moving the jth node from the interpolation, resulting in an interpolation of degree 
n- 1: 

n — 1 

g„-i{x) = ^4''~^^pi{x) 

i=0 
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which still interpolates the integrand at the remaining n nodes. 

The coefficients c^"~^^ of gn-\{x) can be computed, as described in [Gonnet 
2009b], using 



(n-l) 



-1) 



0...n 



where the Cj are the computed coefficients of gn{x) and the are the coefficients 

of the downdated Newton polynomial computed by solving the upper-triangular 
system of equations 



/ ao -{xj +^i) 



V 



72 



an-2 



-{Xj +I3n-l) 
an-1 



-{Xj + I3n) 



OLrx 



I 



/b^\ 
b2 



\bnj 



using back-substitution, where the bi arc the coefficients of the Newton polynomial 
over the nodes of the quadrature rule (Equation (16)) and the ccj, /3j and 7^ are the 
coefficients of the three-term recurrence relation satisfied by the polynomials of the 
orthogonal basis: 

akPk+i{x) = {x + Pk)Pk{x) - 7kPk-i{x). 

The modified vectors c*^""^) and b*^""^) are then used in the same way as c and 
b respectively for the computation of the integral and of the error estimate. 

5. DIVERGENT INTEGRALS 

Divergent integrals are integrals which tend to ±00 and thus cause most algorithms 
to either recurse infinitely or return an incorrect finite result. They are usually 
caught by limiting the recursion depth or the number of function evaluations ar- 
tificially. Both approaches do not per se attempt to detect divergent behavior, 
and may therefore cause the algorithm to fail for complicated yet non-divergent 
integrals. 

Ninham [1966] studied the approximation error when computing 

rh 

x°' dx (20) 



/ 

Jo 



using the trapezoidal rule. The integral exists for a > — 1 and is divergent otherwise. 
Following his analysis, we compute the refined error estimate described in Section 3 
(Equation (18)) for the intervals [0, h] and [0, h/2] using an 11-node Clenshaw-Curtis 
quadrature rulc^" and removing the singular node at x = as described above. 

We note that as the algorithm recurses to the leftmost interval, the local error 
estimate as well as the computed integral itself remain constant for a = —I and 
increase for a < —1. In Figure 2 we plot the ratio of the error estimate in the left 
sub-interval over the error in the entire interval, £[0, h/2]/e[0, h], over the parameter 

^"Remember that we are interpolating over the Chebyshev nodes which is equivalent to using a 
Clenshaw-Curtis quadrature rule. 
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Fig. 2. Ratio of the error estimates for x" dx for the intervals [0, h/2] over [0, h] for different 
a. Note that the error grows (ratio > 1) for a < —1, where the integral is divergent. 

a. For a = —1 the ratio is 1 meaning that the error estimate of the leftmost interval 
remains constant even after halving the interval. For a < — 1, for which the integral 
diverges, the error estimate in the left half-interval is larger than the error estimate 
over the entire interval. 

The rate at which the error decreases (or, in this case, increases) may be a good 
indicator for the convergence or divergence of the integral in Equation (20), where 
the singularity is at the edge of the domain, yet it does not work as well for the 
shifted singularity 

/ |x-/?rdx, /3e[0,h/2]. (21) 
Jo 

Depending on the location of the singularity {x — f3), the ratio of the error estimates 
over [0,h] and [0,h/2] varies widely for both a > — 1 and a < — 1 and can not 
be used to determine whether the integral diverges or not. In Figure 3 we have 
shaded the regions in which the ratio of the error estimates e[0,h/2]/e[0,h] > 1 
for different values of a and the location of the singularity /3. For this ratio to be 
a good indicator for the value of a (and hence the convergence/divergence of the 
integral), the shaded area should at least partially cover the lower half of the plot 
where a < — 1, which it does not. 

A more reliable approach consist of comparing the computed integral in two 
successive intervals [a,b] and [a, {a + b)/2] or [(a + b)/2,b]. For the integrand in 
Equation (20), the integral in the left sub-interval [0, h/2] is larger than that over 
the interval [0, h] when a < —1. For the integral in Equation (21) the ratio of the 
integrals in the intervals [0, h] and [0, h/2] is larger than 1 for most cases where 
a < — 1 (see Figure 4). 

Although this relation (ratio > 1 a < — 1), which is independent of the interval 
h, is not always correct, it can still be used as a statistical hint for the integrand's 
behavior during subdivision. In the area — 2<a<— litis correct approximately 
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Fig. 3. Contour of the ratio of the error estimates for \x — f}]"' dx over the intervals [0, 1] and 
[0, 1/2]. The filled area represent the region in which this ratio is larger than 1. 




Fig. 4. Contour of the ratio of the integral estimates for \x — jSl" dx over the intervals [0, 1] 
and [0, 1/2]. The filled area represent the region in which this ratio is larger than 1. 



two thirds of the time. We wiU therefore count the number of times that 



Qn[a,ia + b)/2] 
Q„[a,b] 



> 1 



,[{a + b)/2,b] 
Qn[a,b] 



> 1 



(22) 



during subdivision. Note that since the integral goes to either +oo or —oo, the 
integral approximations will be of the same sign and thus the sign of the ratios do 
not matter. If this count exceeds some maximum number and is more than half of 
the recursion depth - i.e. the ratio of integrals was larger than one over more than 
half of the subdivisions - then we declare the integral to be divergent and return 
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an error or warning to the user to this effect. 
6. THE ALGORITHM 

We present the aigorithm in two variants, Algorithm 3 and Algorithm 4, using both 
the naive and the refined error estimates presented in Section 3. Both algorithms 
follow the globally adaptive general scheme shown in Algorithm 2. 

The first algorithm (Algorithm 3) uses the naive error estimate in Equation (10) 
in a doubly adaptive strategy using dmax+1 rules of degree n, = 2ni_i, i = 1 . . . dmax, 
using the nodes x^'^ and the Vandermonde-like matrices P('\ i = 0...rfmax- For 
the tests in Section 7, no = 4, hint = 0.1 and dmax = 3 were used. 

In Lines 1 to 4, the coefficients of the two highest-degree rules are computed and 
used to approximate the initial integral and error estimate. In Line 5 the heap 
H is initialized with this interval data. The algorithm then loops until the sum 
of the errors over all the intervals is below the required tolerance (Line 7). At 
the top of the loop, the interval with the largest error is selected (Line 8). If the 
error in this interval is below the numerical accuracy available for the rule used or 
the interval is too small {i.e. the space between the first two or last two nodes is 
zero. Line 11), the interval is dropped and its error and integral are accumulated 
in the excess variables Sxs and gxs (Line 12). If the selected interval has not already 
used the highest-degree rule (Line 14), the coefficients of the next-higher degree 
rule are computed and the integral and error estimate are updated (Lines 15 to 
19). The interval is bisected if either the highest-degree rule has already been 
applied or if when increasing the degree of the rule the coefficients change too 
much (Line 20), analogously to the decision process suggested by Venter and Laurie 
[2002]. For the two new sub-intervals, the coefficients for the lowest-degree rule are 
computed (Line 30) and used to approximate the integral (Line 31). The number 
of times the integral increases over the sub-interval is counted in the variables nr^iv 
(Line 32) and if they exceed nrdivmax and half of the recursion depth nrrec of that 
interval, the algorithm aborts (Line 33) as per Section 4. Note that since the test in 
Equation (19) requires that both estimates be of the same degree, we will use, for 
the estimate q^^^ from the parent interval, the estimate which was computed for the 
first rule. The error estimate for the new interval is computed by transforming the 
interpolation coefficients from the parent interval using Equation (6) and using its 
diflference to the interpolation in the new interval (Line 34) . When the sum of the 
errors falls below the required tolerance, the algorithm returns its approximations 
to the integral and the integration error (Line 40). 

The second algorithm (Algorithm 4) uses the refined error estimate in Equa- 
tion (18), which re-uses the coefficients from a previous level of recursion. For the 
results in Section 7, n = 10 and i?! = 1.1 were used. 

In Lines 1 to 5 an initial estimate is computed and used to initialize the heap 
H. The error estimate is set to oo since it can not be estimated (Line 4). While 
the sum of error estimates is above the required tolerance, the algorithm selects 
the interval with the largest error estimate (Line 8) and removes it from the heap 
(Line 9). As with the previous algorithm, if the error estimate is smaller than 
the numerical precision of the integral or the interval is too small, the interval is 
dropped (Line 10) and its error and integral estimates are stored in the excess 
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Algorithm 3 int_naivc (/, a, 6, t) 

1: for i = . . . do fi^f ((a + b)/2 -{a- b)x["'^ /2j end f 

{evaluate the integrand at : 
2: c(<*-»<-i) ^ (p(rfmax-i))-if[i : 2 : n„ + 1], c('^"'='<) <(- (p('i">ax))-if 



{evaluate the integrand at the nodes x]" ) 

2: c(<*-»<-i) ^ (p(rfmax-i))-if[i : 2 : n„ + 1], c('^"'='<) <(- (p('i">ax))-if 

{com,pute the interpolation coefficients) 

3: qo ^ {b — a)u}^ c^'^'""^'' /2 {approximate the integral as per Equation (5)) 

4: £0 — a)||c(''™'<) — €^'^'""'"■^'11/2 {approximate the error) 

b: H -ir- {[a, 5, c^'''"^"-', go) s^O) c^max) 0) 0]} {init the heap with the first interval) 

6: £xs 0, (/xs ^ {init the excess error and integral) 

7: while X^E.gjyEi > r do 

8; A: ^ argmaxfc 6^ {get the index of the interval with the largest error) 

9: m ^ {uk + bk)/2, h ^ {bk - ak)/2 

10: split false {init split) 

11: if Ek < |5fe|£machCond(P('''=)) V interval too small then 

12: ^ £x5 + £fe, ?xs ^ 9xs + 9fe {collect the excess error and integral) 

13: H -(^ H\ {[ttfc, 6fe, c"'**, gfe, £fe, dk, nrjiv, nrrec]} {remove the kth interval) 

14: else if dk < (imax then 

15: dk ^ dk + 1 {increase the degree in this interval) 

16: for i = . . . ridj, do fi^ f (m- hxf''^/2^ end for 

17: 0^*^*=^ = (p('^fc))~if {compute the new interpolation coefficients) 

18: <s— h(jj^ c^'^''^ {approximate the new integral) 

19: £fe ^ /i||c('^'=) — c°'''|| {approximate the new error) 

20: split ^^'^yjcd^ry > hint {check change in the coefficients) 

21: else 

22: split ^ true (spZii the interval if we are already at highest- degree rule) 

23: end if 

24: if split then 

25: H -(^ H\ {[ok, bk,c°^^, qk,£k, dk, nrjiv) """rec]} {remove the kth interval) 

26: for i = ... no do 

27: ^ / ((afe + m)/2 + /ia;f /2) , Z^'* ^ / ((m + 6fe)/2 + hx^"^ /2) 

28: end for 

29: for half e {left, right} do 

30: c^^"' = (p(o))-if''^'f {compute the new interpolation coefficients) 

31: 9haif ^Cq^'^/v^ {approximate the new integral) 

32: if ghaif > then nr^fjf ^ nrf + 1 else nr^f^f ^ nr°!^ end if 

33: if nr^?)' > nrdivmax A 2nrj?J*^ > nrrec then return error end if 

{abort on divergence) 

34: £haif ^ /illc*^^'^ — T^^'^c"'"^!! {approximate the new error) 

35: end for 

36: H HU {[ak, m, c^^^, gieft, £ieft, 0, nr||jft, nr^ec + 1], 

37: [m, bk, c"^^^, bright, £right, 0, nr^'.f , nrrec + 1]} 

{push the new intervals back on the heap) 

38: end if 

39: end while 

40: return q^ + ^ w </i >. \Sxs + sir £i] {return the integral and the error) 

ACM Transactiojns on JVfatfHegiaticaJ Software, "VMj^NoJ N; Month 20YY. ^ ' 



Adaptive Quadrature Using Explicit Interpolants • 17 



Algorithm 4 int_refined (/, a, b, r) 



for z = . . . n do /i / ((a + 6)/2 - (a - b)xi/2) end for 

{evaluate the integrand at the nodes Xi) 

{compute the interpolation coefficients) 
(7o ^ (6 — a)w^c(''"'")/2 {approximate the integral) 

£o 00 {start with a somewhat pessimistic appreciation of the error) 

H {[a, b, c, b, qo, sq, 0, 0]} {init the heap with the first interval) 

Exs ^ 0, (7x5 ^ {init the excess error and integral) 

while Y^g.^j^Si > r do 

k argmaxfc £fe {get the index of the interval with the largest error) 

H -(^ H\ {[ofe, 6fe, c°''', b"'*^, (/fc, efe. nr^j^. nrrec]} {remove the kth interval) 

if Ek < |(j'/c|emachCond(P) V interval too small then 
^xs ^xs + £fc) 9xs 9xs + Qk {collect the excess error and integral) 

else 

{ak + bk)/2, h ^ {bk - ak)/2 
for i = . . . n do 

fl"^ ^ f {{ak + m)/2 + hx,/2), f^^' ^ f {{m + bk)/2 + hx,/2) 

{evaluate the integrand at the nodes Xi in the sub -intervals) 

end for 

for half € {left, right} do 

^.haif _ (^p^-ifhaif {computc the new interpolation coefficients) 

9haif hu}^c^^^^ {approximate the new integral) 

if ghaif > qk then nr^fjf ^ nr^!^ + 1 else nr^fj' ^ nr^!^ end if 
if nr^?)' > nrdivmax A 2nrj]?j*^ > nr^ec then return error end if 

{abort on divergence) 

fha\t^^ ^ l|bi!i^-2"+iT(i"%b°ii|| {approximate the higher derivative) 

if max||Pc°''' - f'^^'fj - t?i/h!^^^^ |Pb°'''|| > then 

£^haif /illc*^^" — c°'''|| {compute the un-scaled error) 

else 

ffhaif /*/hIif^^'l|t>''^'^|| {compute the extrapolated error) 

end if 
end for 

H^HU{[ak,m, c'^f*, b'^f, gieft, eieft, nrif^, nr.ec + 1], • • • 
[m, bk, c"s^^t, b^'g^S gnght, bright, nr|jt, nr.ec + 1]} 

{push the new intervals back on the heap) 

end if 

end while 

return j^^xs + X^g-g// Qi ) [^xs + X^E^gj^ £i] {return the integral and the error) 



ACM Transactions on Mathematical Software, Vol. V, No. N, Month 20YY. 



18 • Pedro Gonnet 



variables £xs and q^s (Line 11). The algorithm then computes the new coefficients 
for each sub-interval (Line 18), as well as their integral approximation (Line 19). 
If the integral over the s\ib-intcrval is larger than over the previous interval, the 
variable nr^iv is increased (Line 20) and if it exceeds nrjivmax and half of the recursion 
depth nrrec, the algorithm aborts with an error (Line 21). In Line 23 the algorithm 
tests whether the conditions laid out in Equation (19) for the approximation of the 
n + 1st derivative hold. If they do not, the un-scaled error estimate is returned 
(Line 24), otherwise, the scaled estimate is returned (Line 26). Finally, both sub- 
intervals arc returned to the heap (Line 29). Once the required tolerance is met, 
the algorithm returns its approximations to the integral and the integration error 
(Line 32). 

Although the algorithm descriptions in Algorithms 3 and 4 arc quite complete, 
some details have been omitted for simplicity. First of all, when the function values 
are computed (Lines 1, 15 and 27 of Algorithm 3 and Lines 1 and 15 of Algorithm 4), 
it is imdcrstood that previously computed function values at the same nodes, i.e. on 
the edges of the domain for both algorithms or inside the Clenshaw-Curtis rules of 
increasing degree for Algorithm 3, are re-used and not re-evaluated. 

We have also not included the downdatc of the interpolations when NaN or ±lnf 
is encountered. This is done as is shown in Algorithm 5. If, when evaluating 
the integrand, a non-numerical value is encountered, the function values is set to 
zero and the index of the node is stored in nans (Line 4). The coefficients of the 
interpolation are then computed for those function values (Line 6). For each index 
in nans, first the coefficients b of the Newton polynomial over the nodes of the 
quadrature rule are down-dated as per Equation (4) (Line 8). The downdatcd b is 
then in turn used to downdate the interpolation coefficients c as per Equation (4) 
(Line 9). 

Furthermore, to improve memory efficiency, both algorithms maintain at most 
200 intervals in the heap. If this number is exceeded, the interval with the smallest 
error estimate is removed and its integral and error estimates are added to the 
excess variables q^s and £xs respectively. 



Algorithm 5 Interpolation downdate procedure 
1: nans {} (initialize nans) 

2: for i = . . .n do 

3: fi^f{o-+ ^iii(6 — a)) {evaluate the integrand at the nodes Xi) 

4: if fi € {NaN, Inf} then fi ^ 0, nans nans U {i} end if 

{if the result is non-numerical, set the node to zero and remember it) 

5: end for 

6: c P~-^f {compute the initial interpolation coefficients) 

7: for i G nans do 

8: b U~^b {downdate the coefficients of the Newton polynomial) 

9: c c — l^b {downdate the coefficients of the interpolation) 

10; n ^ n — 1 {decrement the degree) 

11: end for 
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7. VALIDATION 

In the following, we will test both algorithms described in Section 6 against the 
following routines: 

— quadl, MATLAB's adaptive quadrature routine [The Mathworks 2003], based on 
Gander and Gautschi's adaptlob [Gander and Gautschi 2001]. This algorithm 
uses a 4-point Gauss-Lobatto quadrature rule and its 7-point Kronrod extension. 

— DQAGS. QlJADPACK's [Pic^ssens et al. 1983] adaptive quadrature routine using a 
10-point Gauss quadrature rule and its 21-point Kronrod extension as well as the 
£-Algorithm to extrapolate the integral and error estimate. The routine is called 
through the GNU Octave [Eaton 2002] package's quad routine. 

— da2glob from Espelid [2007], which uses a doubly-adaptive strategy over rules 
of degree 5, 9, 17 and 27 over 5, 9, 17 and 33 equidistant nodes^^ respectively, 
using the error estimator described in Berntsen and Espelid [1991]. 

The two algorithms in Section 6 were implemented in the MATLAB programming 
language^^. 

Over the years, several authors have specified sets of test functions to evaluate 
the performance and reliability of quadrature routines. In the following, we will 
use, with some minor modifications, the test "families" suggested by Lyness and 
Kaganove [1977] and the "battery" of functions compiled by Gander and Gautschi 
[1998], which are an extension of the set proposed by Kahaner [1971]. 
The function families used for the Lyness-Kaganove test are 

f-i 

|a;-A|"da;, A e [0, 1], a e [-0.5, 0] (23) 



Jo 



J 

Jo 

[ {x> A)e"^ dx, A G [0, 1], a e [0, 1] (24) 

1 

exp(-a|a; - A|)dx, A e [0, 1], a e [0,4] (25) 



(26) 



I 



10"/((a;- A)2 + 10")dx, A e [1,2], a e [-6,-3] 

4 

^10"/((a;-Ai)2 + 10")da;, A» € [1, 2], a S [-5, -3] (27) 

i=l 
L 

2/3(a; - A) cos(/3(a; - A)^) dx, A e [0, 1], a e [1.8, 2], (28) 



/3 = 10°/max{A^(l-A)^} 



where the boolean expressions are evaluated to or 1. The integrals are computed 
to relative precisions'^ of r = 10~^, 10~^, 10~^ and 10~'^ for 1000 realizations of 

^^The final rule of degree 27 over 33 nodes is constructed such as to maximize its numerical 
stability. 

'^^The source-code of both routines is available online at 
http : / /people . Inf . ethz . ch/ gonnetp/ toms/ 

Since all algorithms tested use an absolute error bound, the exact integral times the relative 
tolerance was used. 
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the random parameters A and a. The results of these tests are shown in Table I. 
For each function, the number of correct and incorrect integrations is given with, 
in brackets, the number of cases each where a warning (cither explicit or whenever 
an error estimate larger than the requested tolerance is returned) was issued. We 
consider an integration to be correct only when the returned value is within the 
required tolerance of the exact result. 

The functions used for the "battery" test arc 

h = /o /i4 = C V50e-50--' dx 

/2 = !l{x>Q.2.)&x /i5 = £°25e-25-da; 

h = /o x^'^ da; /le = Jq 50(7r(2500a;2 + l))-i dx 

= /_^(|| cosh(a:) — cos(a;)) da; fn = Jq 50(sin(507ra;)/(507ra:))^ da; 

/5 — ^^^{x'^ + a:^ + 0.9)"-^ da; /is = /J^ cos(cos(x) + 3sin(a;) + 2cos(2a;) + 3cos(3x)) dx 

/e = /Qa;3/2da; /19 = /„Mog(a;) da; 

/7 = Jo da; /20 = £1(1.005 + a;2)-i da; 

h = /o (1 + x^)-^ dx /21 = /o ELi [cosh(20*(x - 2z/10))] da; 

/g = /g 2(2 + sin(107ra;))-Ma; /22 = 47r2a; sin(207r.x) cos(27rx) da; 

/lo = /o (1 + dx /23 = /o (1 + (230a; - 2,Qf)-^ dx 

fn = /o (1 + e")-i da; /24 = /' [e'J da; 

fi2 = Jo x{e^ - 1)"' da; /25 = J^ix + l){x < 1) + (3 - a;)(l < a; < 3) 
/i3 = sin(1007ra;)/(7ra;)da; +2(a;>3)da; 

where the boolean expressions in /2 and /25 evaluate to or 1. The functions are 
taken from Gander and Gautschi [1998] with the following modifications: 

— No special treatment is given to the case a; = in /12, allowing the integrand to 
return NaN. 

— /i3 and fn arc integrated from to 1 as opposed to 0.1 to 1 and 0.01 to 1 
respectively, allowing the integrand to return NaN for x = 0. 

— No special treatment of a; < 10~^^ in /19 allowing the integrand to return —Inf. 

— /24 was suggested by J. Waldvogel as a simple yet tricky test function with 
multiple discontinuities. 

— /25 was introduced in Gander and Gautschi [1998], yet not used in the battery 
test therein. 

The rationale for the modifications of /12, /13, fn and /19 is that we can't, on one 
hand, assume that the user was careful enough or knew enough about his or her 
integrand to remove the non- numerical values, and on the other hand assume that 
he or she would still resort to a general-purpose quadrature routine to integrate it. 
Any general purpose quadrature routine should be robust enough to deal with any 
function, provided by either careful or careless users. 

These changes have little effect on quadl and DQAGS since, as mentioned in Sec- 
tion 4, the former shifts the integration boundaries by £mach if a non-numerical 
value is encountered on the edges of the integration domain and the later uses 
Gauss and Gauss-Kronrod quadrature rules which do not contain the end-points 
and thus avoid the NaN returned at a; = for /12, /13 and fn and the Inf at a; = 
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in /19. da2glob treats NaNs by setting the offending function values to 1. Due 
to the rather fortunate coincidence that /12, /13 and /17 all have a limit of 1 for 
a; — > 0, this integrator is in no way troubled by these integrands. 

The battery functions were integrated for the relative tolerances r = 10^'^, 10^^, 
10~® and 10~^^ and compared to the exact result, computed analytically. The re- 
sults are summarized in Tabic II, where the number of required function evaluations 
for each combination of integrand, integrator and tolerance are given. If integration 
was unsuccessful, the number is stricken through. If the number of evaluations was 
the lowest for the given integrand and tolerance, it is shown in bold face. 

Finally, the integrators were tested on the problem 



for 1000 realizations of the random parameter A and different values of a e 

[— 0.1,-2] for a relative^^ tolerance r = 10""^. Since for a < —1, the integral 
diverges and can not be computed numerically, we are interested in the warnings 
or errors returned by the different quadrature routines. The results are shown in 
Table III and Figure 5. For each integrator we give the number of successes and 
failures as well as, in brackets, the number of times each possible error or warning 
was returned. The different errors, for each integrator, are: 

— quadl: (Min/Max/Inf) 

— Min: Minimum step size reached; singularity possible. 

— Max: Maximum function count exceeded; singularity likely. 

— Inf: Infinite or Not-a-Number function value encountered. 

— DQAGS: (ieri/ier2/ier3/ier4/ier5) 

— ieri: Maximum number of subdivisions allowed has been achieved. 

— ier2: The occurrence of roundoff error was detected, preventing the requested 
tolerance from being achieved. The error may be under-estimated. 

— iers: Extremely bad integrand behavior somewhere in the interval. 

— ier4: The algorithm won't converge due to roundoff error detected in the extrap- 
olation table. It is presumed that the requested tolerance cannot be achieved, 
and that the returned result is the best which can be obtained. 

— iers: The integral is probably divergent or slowly convergent. 

— da2glob: (noise/min/max/sing) 

— noise: The requested tolerance is below the noise level of the problem. Required 

tolerance may not be met. 
— min: Interval too small. Required tolerance may not be met. 
— nnax: Maximum number of function evaluations. Required tolerance may not 

be met. 

— sing: Singularity probably detected. Required tolerance may not be met. 

— Algorithms 3 and 4: (err/div) 

- err: The final error estimate^ is larger than the required tolerance. 
— div: The integral is divergent. 



For a < — 1 an absolute tolerance of r = 10 ^ was used. 




A e [0, 1] 
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-0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6 -1.8 -2 

DQAGS 




-0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6 -1.8 -2 

da2glob 




-0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6 -1.8 -2 

Algorithm 4 




Fig. 5. Results of computing \x — A|™ dx for 1 000 realizations of A £ [0, 1] for different values 
of a (x-axis) . The curves represent the number of correct integrations and the number of times 
each different warning was issued for each value of a. 

Thus, the results for DQAGS at a = —0.8 should be read as the algorithm returning 
146 correct and 854 false (requested tolerance not satisfied) results and having re- 
turned the error \er^ (bad integrand behavior) 58 times and the error iers (probably 
divergent integral) 4 times. 

8. DISCUSSION 

As can be seen from the results in Table I, for the integrands in the Lyness-Kaganove 
test, both new algorithms presented in Section 6 are clearly more reliable than 
quadl and DQAGS. MATLAB's quadl performs best for high precision requirements 
(small tolerances, best results for r = 10~^), yet still fails often without warning. 
QUADPACK's DQAGS does better for low precision requirements (large tolerances, 
best results for t = 10~^), yet also fails often, more often than not with a warning. 

Espelid's da2glob does significantly better, with only a few failures at r = 10"'^ 
(without warning) and a large number of failures for Equation (28) at r = 10"^^, 
albeit all of them with prior warning. The former are due to the error estimate 
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not detecting specific features of ttie integrand due to ttie interpolating polynomial 
looking smooth, when it is, in fact, singular (see Figure 7). The latter were due 
to the integral estimates being affected by noise (most often in the 17-point rule), 
which was, in all cases, detected and warned against by the algorithm. 

The new algorithms fail only for Equation (23) at small tolerances since the in- 
tegral becomes numerically impossible to evaluate (there are not sufficient machine 
numbers near the singularity to properly resolve the integrand) , for which a warn- 
ing is issued. This problem is shared by the other integrators as well. Algorithm 3 
also fails a few times when integrating Equation (27). In all such cases, one of the 
peaks was missed completely by the lower-degree rules, giving the appearance of a 
flat curve. Both algorithms also failed a few times on Equation (28) in cases where 
the resulting integral was several orders of magnitude smaller than the function 
values themselves, making the required tolerance practically un-attainable. 

Whereas the new algorithms out-perform the others in terms of reliability, they 
do so at a cost of a higher number of function evaluations. On average. Algorithm 3 
uses about twice as many function evaluations as da2glob, whereas Algorithm 4 
uses roughly six times as many. 

The results are summarized in Figure 6. The plots, for each tolerance and inte- 
grand, show where each algorithm stand in terms of relative reliability and relative 
efficiency. If we divide the plots into four regions 



iifficicnt and 
reliable 


slow yet 
reliable 


efficient yet 
unreliable 


slow and 
unreliable 



we can sec that whereas da2glob is clearly efficient and reliable, the two new algo- 
rithms are slightly more reliable, yet slow. The results for quadl and QUADPACK's 
DQAGS are scattered over all four regions. 

This trend is also visible in the results of the "battery" test (Tabic II). Algo- 
rithms 3 and 4 fail on /21 for all but the highest and second-highest tolerances 
respectively, since the third peak at a; = 0.6 is missed completely. da2glob also 
docs quite well, failing on /21 at the same tolerances and for the same reasons as 
the new algorithms and on /24 for t < 10~^, using, however, in almost all cases, 
less function evaluations than Algorithms 3 or 4. 

It is interesting to note that quadl, DAQGS and da2glob all failed to integrate 
for tolerances r < 10""^. A closer look at the intervals that caused each algorithm 
to fail (see Figure 8) reveals that in all cases, multiple discontinuities in the same 
interval caused the error estimate to be accidentally small, leading the algorithms to 
erroneously assume convergence. This is not a particularity of the interval chosen: 
if we integrate 

/ [e^Jdx, A €[2.5,3.5] (29) 
Jo 

for 1 000 realizations of the random parameter A for r = 10~^ using these three in- 
tegrators, they fail on 894, 945 and 816 of the cases respectively. Both Algorithms 3 
and 4 succeed in all cases. 
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Fig. 6. Scattcr-plots of the results of the Lyness-Kaganove test-suite for each tolerance. Each 
point represents one of the test functions (Equations 23 to 28). Its location is determined by 
the relative number of function evaluations (on the x-axis) and the relative number of correct 
evaluations (on the j/-axis). 

One could argue that the different error estimates are all equally good and the 
different ratios of reliability vs. efficiency are only due to their parameterization, 
i.e. their scaling of the error estimate. We can test this hypothesis by scaling the 
error estimates of all five algorithms by the smallest value^^ p such that the 1 000 
runs at T — 10"'^ for Equation (23) produce no incorrect results: 



T = 10"^ 


quadl 


DQAGS 


da2glob 


Algorithm 4 


Algorithm 3 




P 


^eval 


P 




P 




P 


^eval 


P 


T^eval 


Eqn (23) 


40 000 


783.18 






50 


229.40 


0.006 


319.24 


0.02 


243.85 



For quadl, a scaling of p = 40 000 was necessary, resulting in an 8- fold increase 
in the number of required evaluations and for QUADPACK's DQAGS, no amount of 
scaling produced the correct result. For da2glob. Algorithm 3 and Algorithm 4, 
only moderate scaling was necessary, bringing the number of required function 
evaluations much closer to each other. There is, therefore, for da2glob a potential 



For simplicity, we consider only values in units of the next-closest power of 10. 
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Fig. 7. Some cases in which the 5, 9 and 17-point rules used in da2glob fail to detect a singularity 
(as in Equation (23)). The assumed integrand (red) is sufficiently smooth such that its higher- 
degree coefficients, from which the error estimate is computed, are near zero. 




Fig. 8. Intervals in which the error estimates of quadl (a), DQAGS (b) and da2glob (c) failed for 
/24. The interpolatory polynomials used to compute the error estimates are shown in red. 
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for tuning empirical scaling factors towards more reliability, yet at the cost of close 
to the same number of function evaluations as Algorithms 3 and 4, at least for 

Equation (23). 

This approach, however, breaks down completely for Equation (29). As can be 
seen in Figure 8, the error estimates which cause the individual algorithms to fail 
are not merely small: they are zero, and hence no amount of scaling will fix them. 
The two new algorithms are therefore not only more reliable due to a stricter (or 
more pessimistic) scaling of the error estimate, but mainly due to the fundamentally 
different type of error estimate, which is less prone to accidentally small estimates 
(see [Gonnet 2009a]). 

In the final test evaluating divergent integrals (Table III), quadl fails to distin- 
guish between divergent and non-divergent integrals, reporting that a non-numerical 
value was encountered for —1.0 > a > —1.4 and then aborting after the maximum 
10000 function evaluations^^ for a < 1.4. For a < —1.0, DQAGS reports the in- 
tegral to be either subject to too much rounding error or divergent. The latter 
correct result was returned in more than half of the cases tested. In most cases 
where a < —0.8, da2glob aborted, reporting that the smallest interval size had 
been reached, or, in some cases, that the maximum number of cvahiations (by 
default 10 000) had been exceeded. All these cases were accompanied by an ad- 
ditional warning that a singularity had probably been detected. For a = —0.8, 
Algorithm 4 fails with a warning that the required tolerance was not met and as of 
a < —1.1 both Algorithms 4 and 3 abort, correctly, after deciding that the integral 
is divergent. 

We conclude that the new Algorithms 4 and 3, presented herein, are more reliable 
than MATLAB's quadl, QUADPACK's DQAGS and Espehd's da2glob. This higher 
reliability is not due to a stricter scaling of the error, but to a new type of error 
estimator which avoids most of the problems observed in the other algorithms. This 
increased reliability comes at a cost of two to six times higher number of function 
evaluations required for complicated integrands such as those in Equations 23 to 
28. 

The tradeoff between reliability and efficiency should, however, be of little concern 
in the context of automatic or general-purpose quadrature routines, which should 
work reliably for any type of integrand. Most modifications which increase efficiency 
usually rely on making certain assumptions on the integrand, e.g. smoothness, 
continuity, non-singularity, monotonically decaying coefficients, etc. . . If, however, 
the user knows enough about his or her integrand as to know that these assumptions 
are indeed valid and therefore that the algorithm will not fail, then he or she knows 
enough about the integrand as to not have to use a general-purpose quadrature 
routine and, if efficiency is crucial, should consider integrating it by a specialized 
routine or even trying to do so analytically. 

In making any assumptions for the user, we would be making two mistakes: 

(1) The increase in efficiency would reward users who, despite knowing enough 

about their integrand to trust the quadrature rule, have not made the eSoit to 
look for a specialized or less general routine. 



'This termination criteria had been disabled for the previous tests. 
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(2) The decrease in reliability punishes users who have turned to a general-purpose 
quadrature routine because they knowingly can not make any assumptions re- 
garding their integrand. 

It is for this reason that wc should have no qualms whatsoever in sacrificing a bit of 
efficiency on some special integrands for much more reliability on tricky integrands 
for which we know, and can therefore assume, nothing. 

Ideally, software packages such as Matlab or libraries such as the Gnu Scientific 
Library [Galassi et al. 2009] should offer both heavy-duty quadrature routines such 
as Algorithms 3 and 4 presented herein, alongside other efficient routines, such as 
da2glob [Espclid 2007] for less trickj^ integrands, as well as even more efficient meth- 
ods {e.g. the exponentially convergent integrator proposed by Waldvogel [2009]) for 
analytic (continuous and smooth) integrands, allowing the user to choose according 
to his or her specific needs or level of understanding regarding the integrand. 
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